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Abstract 

We consider a nonlinear partial differential equation of conservation type to describe the dy- 
namics of vesicular stomatitis virus observed in aliquots of fixed particle number taken from an 
evolving clone at periodic intervals of time Q]. The changes in time behavior of fitness function 
noticed in experimental data are related to a crossover exhibited by the solutions to this equation 
in the transient regime for pulse-like initial conditions. As a consequence, the average replication 
rate of the population is predicted to reach a plateau as a power t~2 . 

PACS: 87.10.-Fe; 02.30.Jr 



1 



The evolution of RNA viruses has been extensively studied in several systems, including 
vesicular stomatitis virus (VSV), foot-and-mouth disease virus (FMDV) and phage — 6 

. This study is facilitated by the fact that virus clones present very high replication and 
mutation rates so a highly heterogeneous population - a quasispecies, as referred in virology 
0] - can be derived from a clone in relative short periods of time. Some of the experiments 
done with VSV populations were directed to study how the process of transmission of fixed 
amounts of pathogens from host to host (much like as in a flu, for example) is expected to 



influence the process of evolution of the fitness of a clone. 

In order to investigate in a theoretical framework some aspects of the interplay between 
these two processes we model the system dynamics, in analogy to the works in Refs. ^ 

n 

and P], by a Markovian process that take place in a properly defined one-dimensional 
(phenotypic) ordered space. The research interest is both conceptual and motivated not 
only by questions in virology but also by their application in different fields, as traffic flow, 
intracellular transport, molecular motors, and others that can be connected to the motion 
of self-driven many-particle systems . We restrict our investigation, however, to systems 
of conservative nature. 

The experimental procedure used in studies of VSV follows the original protocol proposed 



and described in details in Ref. In short, at instant nr of the n — th transmission, 

Ti = 1,2,3,..., a number A^„(0) ~ 10^ of viral particles is taken from an evolving clone (ec) 
of VSV of initial low fitness and inoculated into a recipient containing a known amount of 
non-infected cells. The viruses are allowed to replicate for a period r = 24 hours, reaching 
a total of Nn{T) particles at the end. This accomplishes for one evolutive passage along 
which selection and mutation processes are expected to take place. Then an aliquot with 
A^n+i(0) ~ 10^ particles is taken from the ec and the procedure is repeated using new (non- 
infected) cells for the next passage. The fitness Fn{T) of the ec at the end of the evolutive 
passage n can be expressed as 

F„(r) = = exp <^ / r,{t)dt } = exp {r,{t„)T} (1) 




iVn(O) 

where the growth rate re(t) of the (ec) is assumed a continuous function of time in the 
considered interval so that the second equality follows from the mean value theorem of 
calculus for t„ within the interval nr < tn < {n + l)r. By taking the initial quantity 
in the aliquot Nn{0) = N{0) = N independent of n, one expects to mimic in this way 



the process of infection by transmission of fixed amounts of pathogens. The quahtative 
properties of the population of these individuals at each passage, and in particular the 
value of re(t„), are expected to change considerably from one passage to another. The aim 
of the experimental procedure is to get information on the time behavior of the relative 
quantity r(tn) = Te(tn) — measured with respect to the growth rate r^j of a non-evolving 
wild type (wt) population taken standard. 



121 1 for n = 1,2,... , are 



Two sets of data obtained for r(tn), measured in units of r ^ 

n 

reproduced in Fig. 1 |5] . One of these (stars) shows a fast increasing regime at initial values 
of n followed by a region where r(t„) grows more slowly as n increases. The other set 
(squares) apparently follows just the fast regime at initial times followed by a region where 
data present strong fluctuations. A plateau (not observed in these data) is expected to be 
attained around passage 80, as reported latter in Ref. |6|. Yet, despite the fact that the data 
presents immense fluctuations at large n so the existence (and localization) of a plateau 
could still be debated, one should not expect that the rate of population increasing grows 
indefinitely due to metabolic constraints and limitations of resources j^. 

Figurel 



A few theoretical works have already been devoted to describe the observed dynamics 
of VSV as in Refs. among others. Although the description in does not 

allow for a finite asymptotic limit, it predicts correctly the two-phase regimes, as observed 
in experiments, but conditioned to the presence of random distributions at initial times. 
The expected limit is obtained in Ref. [BI] by a modified version of the model in j^. Our 
motivations to study these data once again from a theoretical perspective came essentially 
from the possibility to investigate a few questions concerning the interplay between the 
dynamics of fitness evolution and that imposed by the particular experimental environment; 
namely (1) what would be an appropriate way to account for effects of adaptability of 
individual components to the particular transmission mode; (2) what are the possible causes 
and in what conditions one should expect two-phase regimes; and (3) how the kinetics of 
fitness is affected by the form of initial population distribution. 

The key point that have guided our formulation is the fact that r(t„) refers to a property of 
the population in the aliquot chosen at the beginning of each passage n. It seems appropriate 
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then to consider r as the unit of time and in connection to this, a process that describes 
fitness evolution but conserves the total number N of individuals. 

To this end, the particles are considered distributed into distinct components (or sub- 
populations) indexed by A; = 0, 1,2, ...each characterized by a fixed fitness value. Let Nn^k 
be the number of particles of the fc-component present in the aliquot at time nr. The 
dynamics of Nn,k from one passage to another is mapped onto a jumping process that take 
place in an appropriate fitness space where to each position k we assign a value kX with 
A a positive constant. Here, kX is defined as the logarithm of fitness of the fc-component. 
Accordingly, kX may depend on its intrinsic growth rate and on the (fixed) time interval 
used to attribute a fitness measure to it. Such dependence, however, does not need to be 
elucidated here. Within this view Nn^k is seen as the occupation number at site k, at time 
nr. The corresponding concentration Cn,k is given by Cn,k = Nn^/N with J2k-^ri,k = N. 
Jumping of particles that can give rise to changes in occupation numbers are assumed to oc- 
cur between neighboring sites at any instant within the interval r between any two passages. 
Moreover, particles at each site are supposed to behave independently from each other; thus 
Nn^k satisfies the following recurrence relation 

Nn+l,k = Nn,k + Pn,k-lNn,k-l + qn,k+lNn,k+l — Pn,kNn,k — qn,kNn,k (2) 

for n = 1,2, .... The quantities Pn,k and qn,k are the transition probabilities between neigh- 
boring sites. Pn,k is the probability for increasing the fitness of the whole population by 
increasing occupation at position [k + 1). q^^k is the probability for decreasing the fitness of 
the whole population by increasing occupation at position {k — 1). The model is properly 
defined when we specify the forms of pn,k and qn,k- We choose 

Pn,k = P - aCn,k (3) 



qn,k = Q + (yCn,k (4) 

Parameters P, Q and a are positive numbers such that 0<P<1, Q = l — P and a < P. 
This choice was inspired by the dynamics of clannish random walkers, proposed by Montroll 
to describe a related process of gas separation In population dynamics, it resembles 
extended logistic models already considered in the literature |3| to study effects 

of competition among the diverse subpopulations. 



Our reasoning is connected to the first question posed above. It comes from observation 

that at each transmission, some of the components of high growth rates may happen to be 

less adapted to the experimental environment so to appear depleted at the time of sampling if 

compared to components of low growth rates. This point, which has already been discussed 

in the literature [18| for considering coinfections and superinfections processes Q, 

emphasizes the necessity to account for the interplay between the intrinsic dynamics of virus 

and that imposed by the particular transmission mode Therefore, in our formulation 

we use the concentrations of the components present in the aliquot at each time (or at the 

time of each passage) as the proper measure of their relative adaptability at that time. The 

strength with which this affects the dynamics of fitness is regulated by parameter a. At the 

genotypic level (microscopic) this parallels effects of random genetic drift by sampling in 

populations of limited size Parameter P is a measure of the fraction of the progeny of 

component with fitness (A; — 1)A present at passage n that in the absence of the terms in a, 

would appear mutated into component of fitness k\ in the following passage. We attribute 

a numerical value to P as characteristic of each clone of VSV, at each state of dilution 

regarding the multiplicity of (cell) infection (m.o.i.). 

We take A — > 0, r ^ keeping the ratio > F finite. In these limits, C„ ^ becomes a 

r 

function C(x, t) of two continuum variables, the position kX ^ x and time nr ^ t for k and 
n ^ oo and Eq. (0) converges to a quasi-linear PDE of conservation type, 

where a = r(2P — 1) and b = AaT are measured in units of t^^. In Eq. (O number 

n 

conservation appears as a local property |20l |. 

The connection of the model with the measurements of the average fitness performed in 
the experiments is made by expressing the exponent in Eq.(IH), up to a constant r„,, as 

r(t„)r = (x) (nr + r) — (x) (nr) (6) 

where in the above limits (x) converges to a continuous function of time. Our proposal 
consists in approaching this by the average 

(x) (t) = / xC{x,t)dx (7) 
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and, consequently 

r{t) = I {x) (t) (8) 
which is the quantity of interest. The solutions to Eq. © for certain initial conditions 

n fi 

C{x,0) can be examined using the method of characteristics |2J|, |22]. We first consider 
a rectangular pulse such that C{x,0) = for xq < A; C{x,0) = h for A < x < A + m 
and C{x, 0) = for x > A + m. The quantities h, m and A are constants, m is the initial 
dispersion and A controls the average(x) (0) at initial time. Normalization of C{x,t) at all 
times is ensured by setting h-m = 1. For this choice of initial conditions, there is occurrence 
of two shocks of characteristics, one at time ts{i) = 0, at Xs(i){0) = A, and another at time 

at position Xs{2){^) = f^i-n 1) + A. This second shock happens when the first shock 

bh 

front encounters the family of characteristics Xr{t) = m + [a — bC] t + A in the region of 
rarefaction that emerges from Xq = A + m. The curves Xs(i){t) = A + [a — b/2m]t and 
Xs(2){t) = at — \j2bt + A + m for t > ts(2) that describe the dynamics of the two shock fronts 
are determined from direct integration of velocities Xs(i), i = 1,2, prescribed by Rankine's 



condition j2^. The solutions C{x,t) = (A + m — x + at)/bt at the rarefaction interpolate 
the regions where C = and C = h. The solution profile C{x,t) for all x is represented in 
Fig. 2 at different instants of time. Observe that because C{x,t) has distinct dependence 
on t, before and after the second shock, the corresponding averages (x) {t) and consequently 
r(t) computed in these different intervals are expected to be distinct functions of time as 
well so as to exhibit a crossover at t = t2s- The mathematical reason for this is, of course, 
the presence of nonlinear terms in a in the evolution equation (jSJ. In fact, before shock, we 
obtain 

r(t) = ri(t) = a - — + — (10) 

reproducing the observed linear behavior of the growth rate of clones of VSV at small times. 
The novel feature comes out from the solutions after shock when r{t) approaches a constant 
equal to a as a power since for t > 
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r{t)=r2{t) = a-^-j= (11) 



We conjecture that such a crossover at the transient regime is in the origin of changes in time 

behavior observed in the evolution of clones of VSV. Notice that for another initial pulse-like 

hx 

(normalized) distribution with a triangular profile, such that C(x, 0) = — for < x < m 

m 

2 

771 

and C{x, 0) = in all other regions, the shock of characteristics occurs at = — - and the 
solutions to r(t) show the same behavior in approaching the constant a at large times. 
Numerical solutions obtained by iterating Eq. (0) for Gaussian initial distributions display 
similar behavior (results not shown). 

Qualitative features exhibited by C{x,t) in Fig. 2 can be regarded as an evidence of a 
compromise between a tendency for maintenance of the structure of the initial distribution, 
that determines the dynamics within the first stage (before shock) and that of spreading 
into many components of relative low concentrations at latter times (after shock). Thus, 
we understand that these results support arguments in the literature according to which 
the first stage of VSV fitness increasing is controlled by mechanisms of natural selection 
acting on the evolving population as a whole whereas the dynamics at the second stage is 
characterized by bottleneck effects acting on low concentration components [6J. Expression 
(Q for t2s suggests that observation of these two regimes, however, is conditioned either 
to the strength of this compromise measured by a (or, b) and to the extent of the initial 
dispersion m. This is in agreement with the fitting of data in Fig.l into ri{t) for one of 
the two sets (square), which present a high value for t2s, and that into r2(t) for the other 
(star) for which results to be relatively small. Thus, the first set appears to be 
localized in the region of linear growth that precedes the shock. The other set (star) seems 
to reproduce the behavior after shock. Different results were discussed in Ref. j8|; according 
to these a condition for observation of two-phase regimes would be the presence of random 
distributions at initial times (not pulse like). 

In conclusion, the model proposed here to describe the interplay between the process 
of transmission and that related to intrinsic fitness evolution allows one to interpret the 
observations in colonies of VSV in terms of a crossover between two time regimes and 
elucidate the conditions for such observation. In addition, the model predicts that the 
asymptotic limit of r{t) is approached as a power t~^/^. The high quality of the fitting 
with r2{t) in Fig.l strongly supports this prediction. This must be contrasted, however, 
with the results in Ref. jsl according to which the temporal behavior of r{t) appears to be 
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determined by a combination of two exponential functions. Hence, a more accurate study of 
this region with emphasis on an interplay between theoretical and experimental views shall 
be crucial to decide on these matters. In turn, this might be relevant to studies in virology 
for predicting novel aspects of the course of evolution of viral populations. 
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Figure Caption 

1 - Experimental data reproduced from Fig. lA MARM clone C (square) and Fig. 2B 
MARM clone D (star) of Ref. P| . The lines represent the best fittings of data using the 
functions in Eq. lfTm) (dashed) and (fTT|) continuum 2^. Inset: the best fitting with r2{t) 
for clone D in approaching the asymptotic regime (few data points were excluded in this 
fitting as an the attempt to minimize fiuctuations and capture the polynomial behavior at 
the relevant region). 

2 - Solution profile for initial rectangular pulse shown at different instants of time, (a) 
At t = 0; (b) at an instant before the second shock (0 < t < t2s)'- C{x,t) = n, for 
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xi{t) < X < X2{t), rarefaction in the region X2{t) < x < X3{t); (c) at an instant after the 

b b 

second shock (t > t2s)- Here xi{t) — A + (a )t, X2{t) — A + m + (a )t and 

2m m 

Xs{t) — A + m + at. 
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Evolutive Passage 
Figure 1 



